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Abstract 

In  the  development  of  more  efficient  and  stable  polymer  electrolyte  fuel  cell  (PEFC),  it  is  important  to  propose  the  optimal  component 
shape  that  can  generate  high  power  and  uniform  the  current  density  distribution  in  a  single  cell.  In  this  study,  our  past  model  was  improved, 
and  simplified  two-dimensional  PEFC  analysis  model  including  flow  and  heat  transfer  of  cooling  water  was  made.  And  PEFC  internal 
phenomenon,  that  is  hardly  measured  experimentally,  could  be  examined  by  using  this  model.  The  influence  of  changing  the  thickness  of 
membrane  and  gas  diffusion  layer  (GDL)  on  the  cell  performance  was  calculated.  As  a  result,  it  was  confirmed  that  it  is  possible  to  improve 
the  cell  output  by  thinning  the  GDL  more  than  the  membrane  in  case  of  low  voltage  and  by  thinning  the  membrane  more  than  the  GDL  in  case 
of  high  voltage,  but  thinning  the  membrane  and  the  gas  diffusion  layer  increased  the  current  density  distribution.  In  addition,  by  arranging 
the  values  of  average  current  density  and  the  current  density  distribution,  the  evaluation  graphs  were  made,  which  became  a  help  of  the  shape 
design  in  the  membrane  and  the  gas  diffusion  layer. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Now,  fuel  cells  are  developed  from  the  viewpoint  of  the 
effective  utilization  of  energy  and  the  protection  of  the  envi¬ 
ronment.  Especially,  polymer  electrolyte  fuel  cell  (PEFC)  is 
expected  as  driving  power  of  vehicles  and  stationary  power 
supply.  And  the  performance  of  PEFC  is  greatly  improved 
because  of  development  of  new  component  and  optimization 
of  the  system.  As  the  PEFC  power  generation  characteristic 
is  affected  by  the  structure,  the  material,  and  the  operating 
condition,  it  is  greatly  important  to  understand  the  correla¬ 
tion  mechanism  of  a  complex  physicochemical  phenomenon 
in  the  cell  in  detail.  But  there  are  hardly  any  researches 
that  examine  PEFC  internal  phenomenon  while  electricity 
is  generated.  As  it  is  very  difficult  to  measure  accurately  the 
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current  density  distribution  and  relative  humidity  distribution 
in  experiments,  the  numerical  analysis  is  an  effective  method 
to  examine  them.  On  the  other  hand,  PEFC  is  demanded  to 
have  long  life  and  high  durability,  and  the  longevity  exam¬ 
ination  is  often  carried  out.  However,  the  accelerated  test 
of  PEFC  is  difficult  because  neither  the  deterioration  factor 
nor  the  mechanism  is  clear.  It  is  thought  that  the  numerical 
analysis  can  be  applied  to  the  investigation  of  such  degrada¬ 
tion  factor  and  mechanism.  Bernardi  and  Verbrugge  [1,2]  and 
Springer  et  al.  [3]  developed  one-dimensional  model  to  the 
direction  of  membrane  thickness,  and  examined  concentra¬ 
tion  distribution  and  water  management  in  PEFC.  Fuller  and 
Newman  [4]  analyzed  developed  two-dimensional  model  to 
the  direction  of  membrane  thickness  and  gas  flow  channel. 
Nguyen  and  White  [5],  and  Yi  and  Nguyen  [6]  developed 
heat  and  water  transport  models  (2-D)  that  accounted  for 
various  operation  condition  and  membrane  hydration  con¬ 
ditions.  On  the  other  hand,  it  is  thought  that  analysis  with  the 
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Nomenclature 

bc  condensation  rate  constant  (Is-1) 

Cj  molar  concentration  of  species  j  (mol  m-3) 

Cp  specific  heat  at  constant  pressure  (J  (kg  K)  - 1 ) 

Eah  the  value  of  reduction  change  of  water  enthalpy 

to  voltage  (V) 

F  Faraday’s  constant  (96485  C  mol-1) 

h  heat  transfer  coefficient  of  gas  (J  (m2  s  K)-1) 

hw  heat  transfer  coefficient  of  cooling  water 
(J  (m2  s  K)-1) 

A//H2o  change  of  water  enthalpy  between  vapor  and 
liquid  (J  mol-1) 

i  current  density  (A  m-2) 

k  thermal  conductivity  of  solid  phase 

(J  (msK)-1) 

ksep  thermal  conductivity  of  separator  (J(msK)-1) 

/d  GDL  thickness  (m) 

/g  gas  channel  depth  (m) 

/w  cooling  water  groove  depth  (m) 

/s  thickness  of  solid  phase  (m) 

/sep  separator  thickness  between  cooling  water  and 
gas  phase  (m) 

Mj  molecular  weight  of  species  j  (kg  mol  - 1 ) 
p  pressure  in  Eq.  (3)  (Pa) 

^H2o,sat  saturated  vapor  pressure  in  stream  (Pa) 
qi  heat  flux  from  solid  phase  to  gas  phase 

(J  (m2  s)-1) 

q2  heat  flux  from  cooling  water  to  gas  phase 
(J  (m2  s)-1) 

<73  heat  value  generated  by  reaction  (J  (m2  s)-1) 

<74  heat  flux  from  gas  phase  to  solid  phase 

(J  (m2  s)-1) 

<75  heat  flux  from  cooling  water  to  solid  phase 
(J  (m2  s)-1) 

<76  latent  heat  value  of  condensation  (J  (m2  s)-1) 

gw(i)  heat  flux  from  gas  phase  to  cooling  water 

(J  (m2  s)-1) 

<7w(2)  heat  flux  from  solid  phase  to  cooling  water 
(J  (m2  s)-1) 

rj  molar  flux  of  species  j  (mol  (m2  s)-1) 

R  gas  constant  (8.314  J  (mol  K)-1) 

7?rea  all  reaction  rate  (Is-1) 

t  time  (s) 

tm  membrane  thickness  (m) 

T  gas  phase  temperature  (K) 

rw  cooling  water  temperature  (K) 

7s  solid  phase  temperature  (K) 

U  overall  heat  transfer  coefficient  between  gas 
and  cooling  water  (J  (m2  s  K)-1) 

If  overall  heat  transfer  coefficient  between  cool¬ 
ing  water  and  solid  phase  (J  (m2  s  K)-1) 
v  flow  velocity  (ms-1) 

V  operation  voltage  (V) 

~x  distance  (m) 


Ax  distance  on  calculation  mesh  to  neighboring 
gas  channel  (m) 

Greek  letters 

a  net  water  transfer  coefficient 

p  density  of  mixture  (kg  m-3) 

Subscripts 
j  species  j 

w  cooling  water 

Superscripts 
a  anode 

c  cathode 

k  anode  or  cathode 

s  solid  phase 

sep  separator 


computational  fluid  dynamics  (CFD)  technique  is  important 
in  order  to  calculate  the  transport  phenomena  in  detail,  and 
such  study  is  increasing  recently.  Um  et  al.  [7]  and  Wang 
et  al.  [8]  have  developed  two-dimensional  model  with  CFD, 
which  included  two-phase  flow.  Dutta  et  al.  [9]  made  a  three- 
dimensional  computational  model  based  on  the  commer¬ 
cial  software  package  Fluent.  Berning  et  al.  [10]  presented 
a  non-isothermal,  three-dimensional  models  and  calculated 
the  distribution  of  current  density  and  concentration  in  the 
straight  channel.  Mazumder  and  Cole  [11]  examined  the  liq¬ 
uid  water  transport  with  three-dimensional  model.  Li  et  al. 
[12]  analyzed  in  small  cell  with  three-dimensional  analysis. 
Um  and  Wang  [13]  compared  the  performance  of  straight 
flow  channel  with  that  of  interdigiteted  flow  channel  by  three- 
dimensional  analysis.  Berning  and  Djilali  [14]  examined  the 
effect  of  the  porosity  and  thickness  of  gas  diffusion  layer  in 
straight  channel  with  three-dimensional  model.  These  PEFC 
numerical  analysis  models  contributed  to  the  optimization  of 
component  design  and  operation  condition,  and  to  the  exam¬ 
ination  of  issues  included  in  present  cell.  However,  these 
studies  calculated  internal  phenomena  in  short  straight  gas 
channel  and  one  serpentine  channel  or  in  small  cell,  and 
there  are  very  few  researches  that  evaluated  various  com¬ 
ponent  shapes  of  actual  size  cell  under  considering  realistic 
calculation  time  and  calculation  resource. 

In  our  past  researches  [15],  the  effects  of  changing  oper¬ 
ation  temperature,  humidify  temperature  and  hydrogen  and 
oxygen  concentration  in  supply  gas  on  the  I-V  character¬ 
istic  of  a  small  PEFC  were  examined  experimentally.  For 
the  experiment,  we  developed  two  models:  one  was  PEFC 
reaction  model  that  could  show  these  influences  on  PEFC 
reaction  characteristics;  the  other  was  PEFC  reaction  and 
flow  analysis  model  that  was  combined  with  the  thermal  flow 
analysis.  With  this  PEFC  reaction  and  flow  analysis  model, 
five  kinds  of  separators  were  evaluated  from  the  viewpoint 
of  gas  flow  condition,  uniformity  of  current  density  and  tern- 
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Fig.  1.  Schematic  diagram  of  PEFC  single  cell. 


perature,  reduction  of  pressure  drop,  and  ejection  of  water.  In 
this  study,  we  made  a  new  simplified  PEFC  analysis  model 
that  included  the  flow  and  heat  transfer  of  cooling  water,  and 
examined  the  influence  of  membrane  and  gas  diffusion  layer 
thickness  on  cell  performance  with  this  model. 

As  the  gas  diffusion  layer  (GDL),  which  is  a  component 
of  PEFC  is  thinner,  the  concentration  overvoltage  is  lower, 
because  the  diffusion  length  becomes  short.  And  the  thinner 
electrolyte  membrane  is,  the  lower  the  resistance  overvoltage 
is,  because  the  ion  transfer  length  becomes  short.  Conse¬ 
quently,  GDL  and  electrolyte  membrane  were  expected  to 
be  thin  in  order  to  reduce  the  cell  voltage  drop.  (Actually, 
those  thicknesses  were  limited  by  the  viewpoint  of  strength, 
manufacture,  etc.)  On  the  other  hand,  from  the  viewpoint  of 
durability  of  catalyst  and  membrane,  it  is  expected  that  the 
current  density  distribution  and  relative  humidity  distribution 
become  uniform  in  overall  cell.  In  this  study  with  numerical 
analysis  method,  the  effect  of  thickness  of  GDL  and  mem¬ 
brane  was  examined  from  the  following  two  view  points;  one 
was  the  achievement  of  high  out  put  density,  and  the  other 
was  the  formation  of  uniform  current  density  distribution  and 
uniform  relative  humidity  distribution. 


2.  Simplified  reaction  and  flow  analysis  of  PEFC 

2.7.  Structure  of  PEFC  analysis  model 

Fig.  1  shows  schematic  of  PEFC  single  cell  that  is  object  in 
this  study.  PEFC  consists  of  a  membrane  electrode  assembly 
(ME A),  two  gas  diffusion  layers  (GDL)  and  two  separators. 
Polymer  membrane  of  MEA  has  proton  conductivity  and 
is  sandwiched  between  two  thin  platinum  electrode  layers. 
MEA  is  sandwiched  between  two  GDL  and  two  separators. 
Table  1  shows  the  specification  of  MEA,  GDL,  and  the  sep¬ 
arator  targeted  in  this  study.  The  flow  channel  shape  of  the 
separator  was  variously  developed.  In  this  study,  we  used  only 
the  serpentine  separator  shown  in  Fig.  2.  The  electrode  area 
was  a  square  with  150  mm  side.  Fifteen  channels,  which  have 


Table  1 

Specification  of  this  PEFC 

MEA 

Thickness  of  membrane 

30,  50,  100  |xm 

Size  of  catalyst  layer 

150  mm  x  150  mm 

Amount  of  Pt 

3.0  gm-2 

GDL 

Size 

150  mm  x  150  mm 

Thickness 

300  pun 

Separator 

Number  of  the  channel 

15 

Channel  width 

1  mm 

Shoulder  width 

1  mm 

Channel  depth 

0.5  mm 

1  mm  groove  width  and  shoulder  width  and  0.5  mm  groove 
depth,  turned  four  times.  Anode  gas  and  cathode  gas  were 
counter  flow.  Cooling  water  flowed  in  the  back  of  the  separa¬ 
tor,  and  the  channel  shape  was  the  same  as  the  gas  separator. 
The  anode  gas  (hydrogen  75%  and  nitrogen  25%  as  reform¬ 
ing  gas)  and  the  cathode  gas  (oxygen  21%  and  nitrogen  79% 
as  air)  were  supplied  to  the  cell,  respectively.  The  supply 
gases  flowed  through  the  humidifier  set  up  on  the  upstream 
of  the  cell,  operating  the  humidifier  temperature  to  control 
the  humidity.  The  inlet  gas  flow  rate  was  set  automatically  by 
hydrogen  utilization,  oxygen  utilization  and  average  current 
density.  The  inlet  cooling  water  flow  rate  was  set  so  that  the 
temperature  of  outlet  cooling  water  can  become  preset  value. 


Fig.  2.  Gas  channel  shape  of  separator. 
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2.2.  Model  development  and  assumption 

First,  the  most  important  point  of  this  model  explanation 
was  shown  as  followings.  The  various  PEFC  numerical  anal¬ 
ysis  models,  which  were  explained  in  Introduction  in  details, 
were  developed  up  to  the  present.  These  analysis  models  were 
summarized  in  followings.  First,  one-dimensional  model  of 
the  direction  across  the  membrane  was  developed.  Next,  two- 
dimensional  model  was  developed  by  adding  the  direction 
along  the  gas  channel  to  one-dimensional  model.  And  then, 
three-dimensional  model  was  developed  by  adding  the  direc¬ 
tion  of  across  the  gas  channel.  On  the  other  hand,  Nguyen  and 
White  [5],  and  Yi  and  Nguyen  [6]  did  not  consider  the  water 
distribution  across  the  membrane  on  the  assumption  that  the 
membrane,  catalyst  layer  and  gas  diffusion  layer  were  too 
thin.  And  they  assumed  that  the  membrane,  catalyst  layer 
and  gas  diffusion  layer  were  unified  homogeneously  and  cal¬ 
culated  the  local  current  density  and  water  transfer  coefficient 
by  using  the  local  concentration  and  the  temperature,  which 
corresponded  to  the  local  current  density  point,  in  both  of 
anode  and  cathode  gas  channels.  The  direction  across  the 
membrane  of  our  model  was  based  on  the  Nguyen’s  model 
and  was  developed  on  the  assumption  that  the  membrane, 
catalyst  layer  and  gas  diffusion  layer  were  unified  homo¬ 
geneously  and  the  local  current  density  and  water  transfer 
coefficient  was  calculated  from  the  local  concentration  and 
the  temperature  in  gas  channel.  The  direction  along  gas  chan¬ 
nel  of  Nguyen’s  model  was  the  plug  flow  model,  and  the 
shape  of  this  model  was  changed  from  the  straight  shape  to 
the  meander  shape  in  this  paper.  Therefore  the  gas  flow  direc¬ 
tion  was  extended  to  the  direction  of  parallel  to  membrane  (v, 
y-direction).  In  the  microscopic  part  including  the  direction 
across  the  membrane,  our  model  may  be  not  as  accurate  as 
other  three-dimensional  models.  But,  the  calculation  speed 
became  faster  and  the  calculation  resources  became  fewer 
than  other  models  by  simplifying.  That  is  why  it  became  pos¬ 
sible  to  calculate  the  influence  on  overall  cell  in  an  actual  size 
cell  with  our  model,  even  though  it  is  difficult  to  calculate  that 
with  other  three-dimensional  detailed  model.  Therefore  our 
model  is  very  effective  to  examine  the  overall  cell  tendency 
in  an  actual  size  cell,  though  the  calculation  of  the  membrane 
thickness  direction  was  sacrificed. 

As  stated  above,  one-dimension  expresses  the  direction 
along  gas  channel  (x  direction),  and  two-dimension  expresses 
the  direction  of  membrane  plane,  which  included  v  direc¬ 
tion  and  the  direction  across  the  gas  channel  (x-y  direction). 
However,  in  the  figure  of  current  density  distribution  and 
relative  humidity  distribution,  x  and  y-axes  were  converted 
to  rectangular  coordinate  in  order  to  help  to  understand 
easily. 

So  far  we  have  been  examining  the  PEFC  reaction  and 
flow  analysis.  In  this  study,  we  improved  the  former  model 
and  developed  the  new  model  that  could  be  used  for  real  cell 
calculation  in  detail. 

Fig.  3  shows  the  PEFC  simulation  model.  As  shown  in  this 
figure,  gas  flow  velocity,  concentration  and  temperature  were 


Back 

plate 


Cooling 
water 

Separator 


Cathode  gas 

Solid  phase 
(MEA&C.DL) 


Cooling 

water 


Temperature:  Tf 


r  Temperature: 

Tc' 

Concentration: 

Cc 

Velocity: 

vc 

,  Pressure: 

Pc\ 

r . 

Current  density: 

N 

/ 

Temperature:  js 

H20  transfer 
coefficient:  a 


Temperature:  Ja 

Concentration:  Ca 

Velocity:  va 

Pressure: _ Pa 

Temperature:  Tf. 


Fig.  3.  Model  for  the  simulation  of  the  PEFC. 


calculated  in  the  gas  channel  on  the  anode  and  cathode.  And  it 
was  assumed  that  the  temperature  distributions  of  MEA  and 
GDL  were  same  as  each  other,  and  that  they  were  unified,  the 
temperature  and  current  density  were  calculated  in  it.  (This 
unified  part  is  expressed  as  the  solid  phase.)  The  governing 
equations  in  this  simulation  were  derived  with  the  following 
assumption. 

1 .  The  gas  flows  uniformly  in  each  channel  of  the  separator. 
Similarly,  the  cooling  water  flows  uniformly. 

2.  The  volume  of  the  condensation  water  is  ignored,  and 
the  water  moves  with  the  gas. 

3.  The  reaction  area  reduction  caused  by  flooding  of  elec¬ 
trode  is  ignored,  and  the  diffusion  prevention  caused  by 
water  condensation  is  ignored. 

4.  Fluid  is  incompressible  Newtonian  fluid  and  ideal  gas. 
Flow  condition  is  laminar  flow. 

5.  Heat  transfer  between  separator  and  gas  is  ignored.  But 
heat  transfer  between  gas  phase,  solid  phase  and  cooling 
water  is  included. 

6.  Cell  voltage  is  uniform  and  constant. 

7.  Only  resistance  overvoltage  and  water  transfer  in  mem¬ 
brane  include  influence  of  temperature. 

8.  In  membrane,  ionic  conductivity,  electro  osmosis  coeffi¬ 
cient  and  water  effective  diffusion  coefficient  that  depend 
on  membrane  humidity  are  determined  by  water  activity 
of  anode  side. 

9.  The  membrane  thickness  and  GDL  thickness  correspond 
to  the  migration  length  of  the  proton  and  water  in  mem¬ 
brane,  and  to  the  diffusion  length  in  GDL,  respectively. 
And  physical  and  electrochemical  properties  of  the  mem¬ 
brane  and  GDL  are  equal  in  any  shape. 

10.  The  gas  crossover  is  disregarded. 
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2.3.  Derivation  of  simplified  PEFC  analytical  model 


In  our  past  method  [15],  in  order  to  calculate  the  gas  flow 
rate  distribution  in  each  channel,  two-dimensional  analysis 
model  was  used.  In  this  study,  One-dimensional  analysis 
as  plug  flow  in  each  channel  was  available  because  of  the 
assumption  that  the  gas  flow  rate  distribution  was  uniform. 
Though  the  separator  shape  was  two-dimension  structure  to 
the  direction  of  the  membrane  face,  the  quasi-two-dimension 
analysis  model  was  made  by  assuming  the  direction  from 
the  inlet  to  the  outlet  to  be  a  positive  v-direction  in  each 
channel  and  by  making  the  channel  meander.  As  a  result,  the 
simplification  of  the  equations  and  the  speed-up  of  the  cal¬ 
culation  became  possible.  However,  in  case  of  calculation  of 
the  solid  phase  temperature  distribution,  it  was  calculated  by 
two-dimensional  analysis  model.  Moreover,  for  the  simpli¬ 
fication  of  the  calculation,  these  three  terms  following  were 
ignored:  the  viscous  term  in  the  motion  equations;  the  heat 
conduction  term  in  the  energy  balance  equations;  the  diffu¬ 
sion  term  in  the  mass  balance  equations. 

The  equation  of  continuity  is  shown  by  the  following 
expression 


dvk 

dx 


where  v  is  the  velocity  of  mixed  gas,  v  is  distance  along  the 
gas  flow  channel,  7?rea  is  all  reaction  rate,  and  the  superscript 
k  is  the  anode  side  or  the  cathode  side.  RYQSL  is  calculated  by 
the  following  equation 


where  /g  is  the  depth  of  gas  channel,  p  is  the  density  of  mixed 
gas,  Mj  is  the  molecular  weight  of  chemical  species  j,  rj  is 
the  reaction  or  condensation  rate  per  unit  area  of  chemical 
species  j. 

The  equation  of  motion  is  shown  by  the  following  expres¬ 
sion 

n?;k 

pk—  =  ~Vpk  +  pkvkRla  (3) 

where  p  is  pressure,  and  the  operator  D/Dt  is  substantial  time 
derivative  that  is  shown  by  the  following  expression 


Dvk  dvk  v  dvk 

- —  - T  v  - 

Dt  dt  dx 


gen,  oxygen,  nitrogen,  vapor  and  condensed  water  in  anode 
and  cathode  channel. 

The  equations  of  energy  are  shown  by  the  following 
expressions 


(Gas)  pkCk 


DTk  q\  +  q\ 


p  Dt 


lk 


+  pkCkTkRkea 


.8^8  ^  t  s vt2 'T’S  ,  <?3  + 


(Solid)  yosCi - =  ksVzr  + 

r  p  dt 


Is 


(Cooling  water)  pwCp(w)  — 


Dt 


lk 

fw 


(6) 

(7) 

(8) 


In  the  energy  equation  of  gas,  Cp  is  specific  heat,  T  is 
temperature,  q\  and  q2  are  heat  fluxes  from  solid  phase  and 
cooling  water,  respectively.  In  the  equation  of  solid  phase, 
k  is  heat  conductivity,  /s  is  the  thickness  of  solid  phase,  q 3 
is  the  heat  value  per  unit  area  the  result  of  electrochemical 
reaction,  ^4  and  q$  are  heat  fluxes  from  gas  and  cooling  water, 
respectively,  q 6  is  latent  heat  flux  of  water  condensation,  and 
the  superscript  s  is  solid  phase.  In  the  equation  of  cooling 
water,  /w  is  the  depth  of  cooling  water  channel,  qw(i)  and 
qw( 2)  are  heat  fluxes  from  gas  and  solid  phase,  respectively, 
and  the  subscript  w  is  cooling  water.  These  heat  flux  and 
heating  value  are  shown  by  the  following  equation 

q\  =  h\r  -  Tk),  q\  =  Uk(Tk  -  Tk), 

qs3  =  (Eah  -  V)i,  q%  =  ha(Ta  -  Ts )  +  hc(Tc  -  Ts), 

qs5  =  Us(a)(Ta  -  Ts )  +  Um(T, l  -  Ts), 

q%  =  -  A//||,o0'|'|2O(1)  +rH20(l))’ 

Al)  =  uk(Tk  -  rw)>  <7w(2)  =  Um(Ts  -  rk)  (9) 


where  h  is  heat  transfer  coefficient  of  anode  or  cathode  gas, 
Eah  is  the  value  of  reduction  change  of  water  enthalpy  to 
voltage,  V  is  voltage,  i  is  current  density,  A//h2o  is  change 
of  water  enthalpy  between  vapor  and  liquid,  Uk  is  overall 
heat  transfer  coefficient  between  anode  gas  or  cathode  gas 
and  cooling  water,  is  overall  heat  transfer  coefficient 
between  anode  or  cathode  side  cooling  water  and  solid  phase. 
These  overall  heat  transfer  coefficients  are  shown  by  follow¬ 
ing  equations 


The  equation  of  chemical  species  j  is  shown  by  the  fol¬ 
lowing  expression 


BE 1 

Dt 


_ j_  _i_  ^k  nk 

/k  ^  j  rea 


(1  /hk)  +  (/seP/kseP)  +  (1  /hky 
£/s(k)  _  _ \ _ 

((/seP  +  /|)/^seP)  +  (l/^) 


(10) 


where  Cj  is  the  concentration  of  chemical  species  j.  The  equa¬ 
tions  of  species  were  derived  to  eight  kinds  of  C^2,  C^2, 
^h20(v)’  ^h20(1  p  ^o2>  C^2,  C^2q(v),  Ch20(1)  that  are  hydro- 


where  /sep  is  separator  thickness  between  cooling  water  and 
gas  phase,  ksep  is  heat  conductivity  of  separator,  h\  is  heat 
transfer  of  cooling  water  of  anode  side  or  cathode  side. 
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The  reaction  and  condensation  rates  of  each  ingredient  are 
shown  by  the  following  equation 


r^2~  2 ~F'  rH20(v)  -  01 J  +  llbc  (  CH20(v) 


pa 

rH20,sat 

RTa 


.a 

N2 


=  0, 


.a 


rH20(l)  —  I  C 


H20(v) 


pa 

/H20,sat 

RT a 


ro2  = 


4  F’ 
+£b. 


H20(v) 


CH20(v) 


=  -(1  +  2a) 


pc 

^H20,sat 

RTC 


IF 


n2 


=  o, 


rH20(l)  -  ~ltbc  (  CH20(v) 


pc 

^H20,sat 

RTC 


(id 


<u 

U 


where  F  is  Faraday’s  constant,  R  is  gas  constant,  FH2o,sat  is 
saturated  vapor  pressure,  a  is  water  transfer  coefficient,  bc  is 
condensation  rate  constant. 

Current  density  i  was  calculated  with  PEFC  reaction 
model  that  authors  built  in  past  study  [15].  The  electromotive 
force,  the  resistance  overvoltage,  the  activation  overvoltage, 
and  the  concentration  overvoltage  were  modeled  from  the¬ 
oretical  equations,  respectively.  On  the  other  hand,  small 
PEFC  with  single  gas  channel  was  made,  which  was  not  much 
affected  by  the  gas  flow  condition  and  the  concentration  dis¬ 
tribution,  and  the  effects  of  changing  the  concentration  of  the 
supply  gas,  the  temperature  of  the  cell,  the  humidifier  tem¬ 
perature,  the  gas  mass  flow,  and  the  gas  channel  length  were 
investigated  by  the  power  generation  examination  of  this  cell. 
By  comparing  the  experimental  results  with  the  calculation 
results  in  all  operation  condition,  the  parameters  of  the  model 
equations  were  decided  with  the  trial  and  error  method.  And 
the  calculation  data  was  able  to  agree  with  the  experimental 
data  by  using  those  parameters.  In  this  model,  the  current  den¬ 
sity  was  calculated  by  the  following  equation  as  the  function 
of  local  concentration  and  temperature. 

i  =  rn  c^2,  caH2 o(v),  cs2,  CCH20(V),  r\  Tc.  r)  (12) 

As  an  example,  when  the  concentration  of  oxygen  in  the 
cathode  supply  gas  was  changed,  the  comparison  between 
the  experimental  data  and  the  calculation  data  of  current 
density-voltage  curve  is  shown  in  Fig.  4.  The  experimen¬ 
tal  data  and  the  calculation  data  agreed  with  each  other  on 
all  conditions.  Similarly,  it  was  confirmed  that  the  I-V  curve 
was  predictable  for  other  conditions  with  this  model. 

It  is  thought  that  water  is  moved  with  two  mechanisms  of 
electro  osmosis  and  back  diffusion  in  the  electrolyte  mem¬ 
brane.  When  one  proton  moves  from  the  anode  side  to  the 
cathode  side,  the  water  movement  coefficient  a  shows  the 
net  number  of  water  molecules  moving  along  with  proton. 
This  is  from  the  method  by  Nguyen  and  White  [5]. 

Eqs.  (1)— (12)  are  connected  in  their  respective  group,  such 
as  reaction,  concentration,  temperature  and  flow,  and  their 
partial  differential  equations  are  discretized  by  finite  differ¬ 
ential  method.  The  boundary  conditions  of  flow  velocity, 


Fig.  4.  Effect  of  02  concentration  on  cell  voltage-current  density  plots  with 
experimental  and  calculated  results. 

temperature,  and  concentration  are  set  as  followings: 

(1)  gas  inlet  boundary:  these  variables  are  constant, 

(2)  gas  outlet  boundary:  the  gradients  of  these  variables  are 
constant. 

Current  density  and  water  transfer  coefficient  were  calcu¬ 
lated  all  over  the  electrode  area.  In  case  of  calculation  of  these 
variables  at  gas  channel  area,  the  local  concentration  and  tem¬ 
perature  at  that  point  were  used  as  the  needed  variables.  In 
case  of  calculation  of  these  variables  at  no  gas  channel  area 
(at  shoulder  area),  the  local  values  at  adjoining  gas  channel 
were  used.  In  this  case,  the  diffusion  length  was  changed 

to  \J (Id)2  +  (Av)2  in  equation  that  calculates  the  limiting 
current  density,  and  gas  diffusibility  to  shoulder  area  was  con¬ 
sidered.  Av  is  the  width  of  calculation  mesh.  And  the  reaction 
rate  of  this  shoulder  area  was  put  in  the  balance  equations  of 
adjoining  mesh  point  at  gas  channel,  and  the  conservation  low 
was  satisfied  in  single  cell.  Those  variables  were  calculated 
until  becoming  stationary  state.  The  relative  errors  of  balance 
equation  of  mass,  species  and  energy  became  1%  or  less  in 
all  the  calculations.  Table  2  shows  the  calculation  system  and 
method  used  in  this  study. 

3.  Results  and  discussions 

In  general,  thinning  membrane  and  GDL  is  effective  in 
the  decrease  of  the  resistance  overvoltage  and  concentration 

Table  2 


Calculation  system  and  method  used  in  this  study 


Analysis  dimension 

One-dimension 

Coordinate  system 

Rectangular  coordinates 

Digitization  of  space 

Finite  difference  method 

Handing  of  convection  term 

QUICKER  method 

Making  time  dispersed 

Euler  scheme 

Flow  analysis  algorithm 

FS  method 

Iteration  of  pressure 

SOR  method 

Mesh 

Square  staggered  mesh 
150  x  150 

14 
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Table  3 


Operation  conditions 


Pressure 

0.1  MPa 

Inlet  gas  temperature 

70  °C 

Inlet  cooling  water  temperature 

70  °C 

Outlet  cooling  water  temperature 

75  °C 

Humidify  temperature 

60,  70  °C 

Operation  voltage 

0.5,  0.6,  0.7  V 

Utilization  rate 

h2 

75% 

o2 

40% 

Inlet  gas  composition 

Anode 

H2,  75.0% 

N2,  25.0% 

Cathode 

02,  21.0% 

N2,  79.0% 

Fig.  5.  Effect  of  thickness  of  membrane  on  current  density  distribution  at 
0.7  V  (average  current  density  (a):  0.276 A  cm-2,  (b):  0.221  A  cm-2,  (c): 
0.146  A  cm-2). 


overvoltage,  respectively.  But,  it  has  not  been  examined  how 
these  influence  on  internal  phenomenon  of  a  real  cell.  Con¬ 
sequently,  in  this  study,  we  examined  the  influence  of  MEA 
and  GDL  thickness  on  current  density  distribution,  relative 
humidity  distribution,  and  cell  performance  with  a  model  that 
can  calculate.  Table  3  shows  operation  condition.  Fig.  5  shows 
the  current  density  distribution  in  case  of  the  membrane  thick¬ 
ness  being  changed  on  0.7  V,  and  Fig.  6  shows  the  anode 
relative  humidity  distribution.  The  solid  line  arrows  show  the 
inlet  and  the  outlet  of  cathode  gas,  the  broken  arrows  show 
the  inlet  and  the  outlet  of  anode  gas  in  Figs.  5  and  6.  In  Fig.  5, 
as  the  membrane  thickness  increased,  the  current  density 
decreased  because  of  the  increase  of  resistance  overvoltage. 
And  the  position  of  the  high  current  density  and  the  low  cur¬ 
rent  density  changed.  In  case  of  30  pan  membrane,  because  of 
the  concentration  overvoltage,  the  current  density  was  high 
in  the  upstream  of  the  anode  and  cathode  channel,  and  it  was 


0.15 
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Fig.  6.  Anode  relative  humidity  distribution  at  0.7  V. 
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Fig.  7.  Relationship  between  membrane  thickness  and  maximum/average 
of  current  density. 


lowing.  In  case  of  thin  membrane,  though  the  current  density 
was  high  on  the  upstream  because  of  the  decrease  in  the  resis¬ 
tance  overvoltage  and  the  maintenance  of  anode  humidity. 
On  the  other  hand,  the  current  density  was  low  on  the  down¬ 
stream  because  of  the  decrease  in  the  oxygen  concentration 
that  results  from  the  consumption  for  reaction.  Moreover,  as 
the  voltage  lowers,  the  influence  of  the  membrane  thickness 
on  the  current  density  distribution  became  small  because  the 
concentration  overvoltage  increased.  Because  the  increase 
of  the  current  density  distribution  influences  on  local  degra¬ 
dation  and  the  durability  of  the  membrane,  it  is  necessary 
to  decrease  this  distribution.  Accordingly,  it  is  important  to 
design  the  cell  having  high  power  density  and  uniform  cur¬ 
rent  density  distribution  as  much  as  possible.  Fig.  8  shows 
the  influence  of  the  membrane  thickness  on  the  ratio  between 
the  cathode  relative  humidity  and  the  anode  relative  humidity 
(cathode  relative  humidity/anode  relative  humidity).  By  the 


low  in  the  downstream.  In  case  of  50  pirn  membrane,  because 
the  influence  of  the  resistance  overvoltage  increased,  the  cur¬ 
rent  density  was  high  at  the  position  of  the  inlet  where  the 
anode  vapor  concentration  was  high.  In  the  place  where  the 
anode  gas  flowed  from  inlet  a  little,  the  anode  vapor  con¬ 
centration  decreased  by  the  electro-osmosis  effect,  and  the 
current  density  decreased.  However,  at  the  position  from  the 
middle  reached  to  the  downstream,  because  of  the  back  dif¬ 
fusion  from  the  cathode,  the  anode  vapor  concentration  was 
high,  and  the  current  density  rose  gradually.  This  tendency 
was  more  remarkable  in  case  of  100  jam  membrane.  It  was 
possible  to  confirm  these  things  from  the  anode  side  relative 
humidity  distribution  of  Figs.  6  and  7  shows  the  influence 
of  the  membrane  thickness  on  the  ratio  between  the  maxi¬ 
mum  current  density  and  the  average  current  density  (maxi¬ 
mum  current  density/average  current  density).  The  rate  was 
increasing  by  thinning  the  membrane,  and  this  caused  like  fol- 


Membrane  thickness  [pm] 


Fig.  8.  Relationship  between  membrane  thickness  and  the  relative  humidity 
ratio  of  anode  and  cathode. 
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Fig.  9.  Effect  of  thickness  of  GDL  on  current  density  distribution  at 
0.5  V  (average  current  density  (a):  0.434  A  cm-2,  (b):  0.415  A  cm-2,  (c): 
0.305  A  cm-2). 


16 


G.  Inoue  et  al.  /  Journal  of  Power  Sources  154  (2006)  8-17 


C/3 

C/3 

0) 

c 

o 


100 


q  80 


<D 

C 

i_ 

Xi 

e 

r- 

s 


(a) 


E 

C/3 

C/3 

0) 

c 

o 


0) 

c 

n 

•_ 

X) 

E 

<D 


60 

40 

20 

10 


100 

80 

60 

40 

20 

10 


10 


(C) 


200  400  600  800 

GDL  thickness  [jam] 


;  ()()() 


100 


I  80 


C/3 
C/3 

O 

•I  60 


<D 

C 

£3 

•— 

X3 

r— 

5 


40  r 


1.25 


1.6 

1.55 

>»  _L 

1 .5  £  ^ 

c  — 

O  >3 

,  "3  C 

1.45  -  y 

C  ^3 

1.4  I  | 

O  C 

E  3 

3  V 
p  Zi. 


1.35 

1.3 

1.25 

1.2 


*  5 

J  > 
^  * 


(b) 


200  400  600  800 

GDL  thickness  [jam] 


1000 


100 


(d) 


200  400  600  800 

GDL  thickness  [jam] 


1000 


Fig.  10.  Effects  of  membrane  and  GDL  thickness  on  average  current  density  and  current  density  distribution  at  0.5  V  and  0.7  V:  (a)  average  current  density  at 
0.7  V,  (b)  current  density  distribution  at  0.7  V  (c)  average  current  density  at  0.5  V,  (d)  current  density  distribution  at  0.5  V. 


drop  of  voltage,  this  ratio  increased,  because  the  amount  of  the 
water  generation  and  electro-osmosis  increased.  Moreover, 
in  case  of  thick  membrane,  this  ratio  increases  because  the 
amount  of  back  diffusion  of  water  decreased. 

Next,  the  influence  of  the  GDL  thickness  was  examined  by 
the  calculation.  Fig.  9  shows  the  current  density  distribution  in 
case  GDL  thickness  is  changed  at  0.5  V.  The  solid  line  arrows 
show  the  inlet  and  the  outlet  of  cathode  gas,  the  broken  arrows 
show  the  inlet  and  the  outlet  of  anode  gas  in  this  figure.  The 
concentration  overvoltage  decreased  by  thinning  the  GDL 
thickness,  and  the  current  density  increases  all  over  the  cell. 
Moreover  in  the  low  voltage,  the  current  density  decreased 
along  the  cathode  gas  flow  pattern  because  the  influence  of 
oxygen  concentration  overvoltage  was  strong. 

Next,  the  above  calculation  result  of  changing  the  mem¬ 
brane  thickness  and  the  GDL  thickness  was  arranged.  Fig.  10 
is  the  evaluation  graph  of  the  membrane  thickness  and  the 
GDL  thickness  concerning  average  current  density  and  cur¬ 
rent  density  distribution.  Comparing  Fig.  10(a),  which  is 
average  current  density  at  0.7  V,  with  (c),  which  is  average 
current  density  at  0.5  V,  it  was  confirmed  that  the  cell  out¬ 
put  could  be  improved  by  thinning  the  GDL  more  than  the 
membrane  in  case  of  low  voltage  and  by  thinning  the  mem¬ 
brane  more  than  the  GDL  in  case  of  high  voltage.  Moreover, 
these  evaluation  graphs  could  be  applied  to  design  the  PEFC 
components.  For  example,  in  case  of  development  of  the  sat¬ 
isfying  cell  the  requirement,  that  the  average  current  density 


is  larger  than  a  certain  value,  and  that  the  current  density 
distribution  is  smaller  than  a  certain  value,  with  these  eval¬ 
uation  graphs,  it  can  be  decided  that  the  size  of  membrane 
and  GDL  in  intersection  of  two  areas  where  those  conditions 
are  satisfied  is  the  best.  Fig.  10  is  the  result  only  about  the 
influence  of  the  shape  of  the  electrolyte  membrane  and  the 
GDL  to  the  average  current  density  and  current  density  dis¬ 
tribution.  As  there  are  many  kinds  of  the  operating  condition 
of  PEFC,  (for  example  cell  and  gas  temperature,  humidifica¬ 
tion  condition,  component  of  inlet  gas,  hydrogen  and  oxygen 
utilization,  operating  voltage,  and  pressure),  it  is  difficult  to 
examine  in  all  these  operation  condition  and  to  arrange  these 
results.  The  specific  best  shape  of  membrane  and  GDL  could 
not  be  described  still  in  this  study.  However,  it  was  thought 
that  the  course  for  the  best  design  of  the  PEFC  component 
could  be  obtained  by  calculating  in  various  conditions  by 
this  PEFC  numerical  analysis  model,  and  by  making  such 
an  evaluation  graph  according  to  the  demand  condition  of 
each  use  and  operation  condition.  On  the  other  hand,  it  will 
be  necessary  to  search  for  the  factor  that  controls  the  output 
performance  and  the  durability  of  the  cell. 

4.  Conclusion 

We  have  improved  our  past  model,  and  presented  the 
more  realistic  PEFC  analysis  model  in  this  study.  The  PEFC 
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phenomenon  was  hardly  measured  experimentally,  but  our 
improved  model  can  examine  it  with  the  numerical  analysis 
in  detail.  The  simulation  shows  that  the  cell  output  and  the 
current  density  distribution  increase  by  thinning  the  mem¬ 
brane  and  the  GDL.  In  addition,  the  evaluation  graph  that 
became  a  help  of  those  shape  designs  was  made  concerning 
the  average  current  density  and  the  current  density  distribu¬ 
tion. 
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